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Abstract 

We continue the study undertaken in Efroimsky (2005a) where we explored the in- 
fluence of spin- axis variations of an oblate planet on satellite orbits. Near-equatorial 
satellites had long been believed to keep up with the oblate primary's equator in the 
cause of its spin-axis variations. As demonstrated by Efroimsky & Goldreich (2004), 
this opinion had stemmed from an inexact interpretation of a correct result by Gol- 
dreich (1965). Though Goldreich (1965) mentioned that his result (preservation of 
the initial inclination, up to small oscillations about the moving equatorial plane) 
was obtained for non- osculating inclination, his admonition has been persistently 
ignored for forty years. 

It was explained in Efroimsky & Goldreich (2004) that the equator precession 
influences the osculating inclination of a satellite orbit already in the first order over 
the perturbation caused by a transition from an inertial to an equatorial coordinate 
system. It was later shown in Efroimsky (2005a) that the secular part of the inclina- 
tion is affected only in the second order. This fact, anticipated by Goldreich (1965), 
remains valid for a constant rate of the precession. It turns out that non-uniform 
variations of the planetary spin state generate changes in the osculating elements, 
that are linear in \ ft \ (where ft is the planetary equator's total precession rate, 
rate that includes the equinoctial precession, nutation, the Chandler wobble, and 
the polar wander). 

We work out a formalism which will help us to determine if these factors cause 
a drift of a satellite orbit away from the evolving planetary equator. 



1 The scope of the project 
1.1 The motivation 

Calculation of the obliquity of a planet (Ward 1973; Laskar Sz Robutel 1993; Touma & 
Wisdom 1994) are always obtained within a simplified model based on representation of the 

* By "precession," in its most general sense, we mean any change of the direction of the spin axis of the 
planet - from its long-term variations down to nutations down to the Chandler wobble and polar wander. 
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planet by a symmetrical rigid rotator, with no internal structure or dissipative phenomena 
taken into consideration. This model yields, via the Colombo (1966) equation, the history 
of the planet axis' inclination in an inertial frame. Thence the evolution of the obliquity 
can be found. The Colombo (1966) equation was derived for a rigid planet in the principal 
rotation state. These assumptions raise questions when it comes to real physics. First, 
a planet is deformable and, thereby, is subject to solar tides. It also tends to yield its 
shape to the instantaneous axis of rotation. (This phenomenon is always acknowledged 
in regard to the Chandler wobble, but never in regard to the equinoctial precession.) 
Second, a forced rotator is never in a principal spin state, and its angular-velocity vector 
is always slightly off its angular-momentum vector. These three phenomena influence the 
equinoctial precession and, through it, the obliquity variations. On the one hand, these 
phenomena are feeble; on the other hand, we are interested in their accumulation over 
the longest time scales, and therefore we are unsure of the outcome. Last, and by no 
means least, the Colombo description of the equinoctial precession ignores the possibility 
of planetary catastrophes that might have altered the planet's spin mode. 

It would be good to develop a model-independent check of whether the planet could 
have maintained, through its entire past, the same equinoctial precession as it has today. 
Such a check is offered by Mars' two satellites. The present proximity of both moons to the 
Martian equatorial plane is hardly a mere coincidence. Hence, the question becomes: could 
Mars have maintained equinoctial precession, predicted by the Colombo model, through 
its entire history without pushing an initially near-equatorial satellite too far away from 
the equatorial plane? 

1.2 The objective 

If it turns out that the equinoctial precession, predicted by the Colombo (1966) model, 
does not drive the satellites away from the equator, or drives them away at a very slow 
rate, then this will become an independent confirmation of this model's applicability. If, 
however, it turns out that the predicted precession of the spin axis leads to considerable 
variations in the satellite inclination relative to the equator of date, this will mean that 
the Colombo model should be further improved or/and that a planetary catastrophe may 
have altered Mars' spin state. 

According to Goldreich (1965) and Kinoshita (1993) the inclination of a near-equatorial 
satellite only oscillates about its initial value, provided the equinoctial precession is uni- 
form. However, even within the simple Colombo model, the equinoctial precession is 
variable. Besides, in these works non-osculating elements were used, circumstance noticed 
by Goldreich (1965) but missed by many authors who employed and furthered his result. 

Whenever the disturbance depends upon velocities (like a transition from inertial axes 
to ones co-precessing with the planet), a mere amendment of the disturbing function 
makes the planetary equations render not the osculating but the so-called contact orbital 
elements whose physical interpretation is nontrivial (Efroimsky & Goldreich 2004). To 
furnish osculating elements, the equations should be enriched with extra terms, some of 
which will not be additions to the disturbing function. 1 Some of them will be of the first 

1 These terms will complicate both the Lagrange- and Delaunay-type equations. The Delaunay equa- 
tions will no longer be Hamiltonian. This parallels a predicament with the Andoyer elements used in 
the theory of rigid-body rotation with angular-velocity-dependent perturbations (Efroimsky 2007; Gurfil, 
Elipe, Tangren, & Efroimsky 2007) 
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order in the velocity-dependent perturbation, others of the second. For uniform precession, 
the first-order extra terms average out, except for a term showing up in the equation for 
dM /dt (Efroimsky 2005a), as predicted by Goldreich forty years ago. Thus, if we address 
not the elements per se but their secular parts, Goldreich's result obtained for the contact 
elements stays also for the osculating ones: the orbit will oscillate about the uniformly 
moving equator but will not shift away from it. 

As demonstrated in Efroimsky (2005a), under variable precession of the spin axis, the 
secular parts of these precession-caused first-order terms are of the first order. Accordingly, 
the secular parts of the osculating elements may differ from those of the contact ones 
already in this order. 

To understand if Mars could have kept through its entire past the same equinoctial 
precession, we need to determine if the satellite orbits might have shifted away from 
the equator in the cause of nonuniform precession. To see how the secular parts of the 
osculating elements evolve, we shall build the required mathematical formalism based on 
the averaged equations. 



1.3 The means 

The motion of a satellite about a precessing oblate planet is most naturally described 
with orbital elements defined in a co-precessing (equatorial) coordinate system. It is also 
convenient to choose the elements to be osculating. The physical interpretation of such 
orbital variables will be most straightforward. 



1.3.1 Exact planetary equations 

The above defined setting is the two-body problem disturbed by two perturbations - the 
gravitational pull of the equatorial bulge and the transition to a non-inertial frame of 
reference associated with the precessing planetary equator. Together, they generate the 
following variation of the Hamiltonian (Efroimsky & Goldreich 2003, Efroimsky 2005a,b): 

AH W = - [ RoUateH + M • (/ X g) + (fl X f) ■ (f2 X f) ] , (1) 

where the oblateness-caused disturbing potential is 



RoUateM = ^ [ 1 " 3 Sin2 % Sill V + ") ] » ( 2 ) 



p e being the mean equatorial radius of the planet, and v standing for the true anomaly. 
The vector 

^ A A A A din A dh/<p T A dh"n T / 1-» \ 

H = fi x x + fi2 y + a*3 z = x + y sm I p + z cos I p , (3) 

is the precession rate of the planetary spin axis. (In the astronomical literature it is 
sometimes referred to as the rotational vector of the equator.) Angles I p and h p are the 
inclination and the longitude of the node of the planetary equator of date relative to that 
of epoch; unit vectors x , y , z denote a coordinate system fixed on the moving equatorial 
plane of date, z being orthogonal to the equator-of-date plane, and x pointing toward 
the ascending node of the equator of date relative to the one of epoch. For details of 
calculation of I p and h p see subsection 2.2.2 and the Appendix. 



3 



Notations / and g stand for two auxiliary vector functions which play an important 
role in the theory. These are the implicit functional dependencies of the unperturbed 
(two-body) position and velocity upon time and six orbital elements. These dependencies 
emerge as a solution 



r = /(d,..,^,*) 

v = gfCx,..,^*) , 
to the reduced two-body problem 



(4) 

df 

dt 



» Gm r 

r + —r Z = (5) 



and define, geometrically, a Keplerian ellipse or hyperbola parameterised with some set of 
six independent orbital elements C, which are constants in the absence of disturbances. 
Under perturbation, these elements are endowed with time dependence. 

This way, our Hamiltonian variation A7i <osc * ) , too, becomes, through composition, 
a function of time and the same six orbital elements used in (4) (these could be the 
Keplerian or Delaunay or Poincare or Jacobi elements). The Hamiltonian variation is 
equipped with superscript "(osc)" in order to emphasise that this is the form taken 
by the Hamiltonian expressed as a function of osculating orbital elements. This clause, 
seemingly trivial and therefore redundant, turns out to be crucially important. As pointed 
out in Efroimsky & Goldreich (2004) and explained in great detail in Efroimsky (2005a), a 
naive development of the planetary equations in precessing frames leads to a Hamiltonian 
variation different from (1); but that Hamiltonian variation tacitly turns out to be a 
function of non-osculating orbital elements. This tacit loss of osculation in problems with 
velocity-dependent perturbations is an old pitfall in orbit calculations. Though some 40 
years ago Goldreich (1965) warned of these difficulties, the issue has until recently been 
ignored in the literature. This has led to a whole sequence of erroneous results. (As 
we already mentioned above, a very similar situation emerges in the rigid-body attitude 
dynamics.) 

For some general-type parameterisation of the instantaneous conies through six orbital 
variables C\ , . . . , C% , the variation-of-parameters equations will read 



dC n 




(6) 



provided these conies are chosen to be always tangent to the trajectory, i.e., provided the 
parameters are chosen to be osculating 2 (Efroimsky & Goldreich 2004, Efroimsky 2005a). 



2 Had we simply amended the Hamiltonian by the above variation AH^ osc ^ , without inserting the 
extra /2-dependent terms into (6), such equations would yield non-osculating elements, ones parametrising 
a family of non-tangent conies. This would happen because the Hamiltonian perturbation depends not 
only upon positions but also upon the canonical momenta. Another way of getting into this hidden trap 
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A more convenient representation of the above equation will be achieved if one includes 
the - (p, x f) d(jlxf)/dC n term in the Hamiltonian: 



[C n Ci 



da 



d "AW 



+ £ 



df ? dg 

x g - / x 



dt dC n ^ \dC n ° J dC n 

the amended "Hamiltonian" being defined through 

"AW = - 



- AH/* 



df_ 

dC n 



(7) 



1 ^ 



Roblate(v) +/*•(/ X g) + -(p,Xf)-(jlxf) 



(8) 



Here the quotation marks are necessary to emphasise that "AW is not the real Hamil- 
tonian variation but merely a convenient mathematical entity. Under this convention, 
and under the assumption that the parameterisation is implemented through the Kepler 
elements, (7) yields the following system of Lagrange-type planetary equations: 



da 



2 

na 



d ( - "AW ) 

dMn 



A* • / x 



df 



(9) 



de 1 — e 2 



dt na 2 e 



d ( - "AW 

dM 



(10) 



(1 - e 2 ) 1 ' 2 
n a 2 e 



doj 
~dt 



cos 2 



na 2 (l — e 2 ) 1 / 2 sin 



(9 ( - "AH") + /9/ xg _ / - x ^g] _~(f x ?£ 



(11) 



+ 



(1-e 2 ) 1 / 2 
na 2 e 



9(-"AW) , /a/ , - 

— ai — + & xg " /x 



<9e 



is to start with the Cartesian or spherical coordinates and momenta, and to perform the Hamilton-Jacobi 
operation. The resulting variables Cj will then come out canonical and will be the well-known Delaunay 
elements. In case the Hamiltonian perturbation depends upon the momenta, these Delaunay elements will 
be non-osculating, i.e., will parameterise a sequence of instantaneous conies non-tangent to the physical 
orbit. (Efroimsky & Goldreich 2003; Efroimsky 2005a) 



5 



di 
dt 



cos % 



na? (1 — e 2 ) 1 / 2 sin % 



d ( - "AH" ) + _ (df x - x 9 
<9cj ^ I <9cj ^ <9w 



- A*- / x 



9/ 

<9cj 



(12) 



na 2 (1 — e 2 ) 1 / 2 sin z 



d(-"AH") _ df _ r dg\ u r df 



na 2 (1 — e 2 ) 1 / 2 sin z 



9 ( - "AH" ) 

9i 



dM 



na? e 



di 



de ^ \ de ^ <9e I ^ I <9e 



(14) 



2 

na 



g(-"A«") ^ /a/ . „- 

Fa + M ' * Xg - /X 



<9a 



A*- / x 



9/ 

<9a 



where terms p, ■ ( (df/dM Q ) x g — (dg/dM a ) x / ) have been omitted in (9 - 10), be- 
cause these terms vanish identically. (For technical details see the Appendix to Efroimsky 
(2005a).) 

1.3.2 The Approximation 

To obtain the first-order (over ft ) secular parts of the osculating elements, we shall carry 
out two operations: 

1. First, we shall throw out 0(ft 2 ) contribution to "AH" and shall assume that 
preservation of the first-order terms and neglect of the second-order ones in the equations 
makes these equations render solutions valid in the first order. This assumption should re- 
main valid for some interval of time, an interval whose actual duration can be determined 
only through accurate numerical simulation. In our analytical developments we shall hope 
that this interval is sufficiently long. 3 



3 In (9 - 14), the quadratic in ft terms in the right-hand side will be of order \ft\ 2 /n. According to 
Ward (1973), the range of values of \p\ for Mars hardly ever exceeded 10~ 3 yr^ 1 . The value of n , for 
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2. Second, we shall substitute both the disturbing function ( — "A7i" ) and the other 
precession-generated (i.e., /x-dependent) terms with their orbital averages. To be more 
exact, the precession rate ft and each of the elements will be considered as a function of 
the true anomaly v , and will be expanded into a Fourier integral which will then be split 
into two pieces - an integral over the band of frequencies less than the orbital frequency 
and an integral over the higher frequencies: 

C 3 = (C,-) + Cf , ft = (ft) + ft" (15) 

The first piece ( (ft) or (C,) ) will be regarded as the secular part, while the second piece 
( ft or Cj ) will be averaged out. One of the outcomes of this treatment will be that the 
left-hand side of the averaged planetary equations will now contain not the time derivatives 
of the elements but the derivatives of their secular parts. To understand the structure of 
the averaged right-hand sides, let us consider some product A(v) B(v) , where A and 
B denote some of the orbital elements or the projection fi± of the vector ft onto the 
instantaneous normal to the satellite orbit (this projection will appear many times in our 
analysis of the planetary equations): 

AB = ((A) + A Q ) ((B) + B*) = (AB) + (A Bf . (16) 
The secular and high-frequency components of this product will read, correspondingly, as 
(AB) = (A) (B) + (A*B*) (17) 

and 

(ABf ee (A) B^ + (B)A* + (jPB* - (A* B*)) . (18) 

An obvious circumstance is that the secular part of the product consists not only of the 
product of the secular parts of the multipliers but also of the term (A^ B^) containing 
resonances. A less evident but crucially important circumstance is that, technically, the 
above separation of timescales is never implemented exactly (unless one deals from the 
very beginning with the Fourier expansions of all the functions involved). Therefore, the 
(imperfectly calculated) high-frequency parts A^ , B^ , and {A B)^ are unavoidably 
contaminated with the lower-frequency modes, modes whose effect may considerably ac- 
cumulate at large times and exert "back-reaction" upon the secular part of the product. 
(Laskar 1990) 

1.3.3 The planetary equations for the first-order secular parts of the osculat- 
ing elements 

Naively, the afore proposed approximation will lead us to a new system of planetary equa- 
tions. It will be identical to the system (9 - 14), except that now the letters a, e, u , f2 , i, M Q 

the Martian satellites, is of order one day -1 . If we now look, for example, at (12), we shall see that the 
quadratic in fl terms surely cannot contribute to di/dt more than an angular degree over a million of 
years, and are quite likely to remain insignificant over dozens of millions years. Whether these terms may 
be neglected at timescales longer than 100 millions of years - should be checked through a comparison of 
our semianalytical model with a comprehensive numerical simulation. As demonstrated in Gurfil, Lainey, 
& Efroimsky (2007), the answer to this question turns out to be positive: the model remains surprisingly 
exact for, at least, 20 Myr. 
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will denote not the osculating elements but their secular parts. Similarly, p, will now 
stand for the secular part of the precession rate. The Hamiltonian will now be substituted 
with 



AH (6//) = - \{R oh late) + (/2'(/xg) 



(19) 



GmJ 2 p e 3 cos 2 i — 1 

4 ^ (1 _ e 2)3/ 2 



Gma (1 — e 2 ) (/ii sini sinf2 — /i 2 sini cosf2 + /j 3 cos 2) 



where, once again, all letters denote not the appropriate variables but their averages. 

By doing so, we would, however, ignore that the /7-dependent terms in (9 - 14) 
contain products of high-frequency quantities (such as, for example, the product of the 

— * — * 

true-anomaly-dependent expression (df/du) xg - (dg/du) x / by the high-frequency 
part of ft in formula (10) ). Averages of such products will contribute to the secular parts 
of the right-hand sides of the approximate planetary equations, as in the example (17). 
(As we shall see below, these inputs will be due to the commensurabilities between the 
orbital motion of the satellite and the short-term nutations of the primary.) Keeping this 
in mind, we should approximate the exact planetary equations rather with the following 
system: 



da 



2 

na 



-<P If 



x 



df 



) 



(20) 



de 
dt 



no? e 



df 



(21) 



3 2U/2 



na? e 



duo — cos i 

dt na?{\ — e 2 ) 1 / 2 sin i 



di \ di di 



+ 



'1 - e 2 ) 1 / 2 



na? e 



(22) 



de 



~, , df _ - d a 
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(23) 



na 2 (1 — e 2 ) 1 / 2 sin i 



d ( - AH( ef V 



on 



+ (p 



df 



x 



<9g 



- ( /* / 



x 




dfi _ 1 

na 2 (1 — e 2 ) 1 / 2 sin i 



d(-AH^) (df 
di \ oi 




{'H [fx % )> , (24) 



dM 



na 2 e 



d(-AH^) (df 
de ^ ^ ^ \ de 




<*f/xg|) 



(25) 



2 

na 



fl(-Att«/>) / 0/ 
da + ( M Ua 




the angular brackets denoting the secular parts. In (16), (17) and (19) we took into account 
the fact that the averaged and truncated Hamiltonian (15) depends neither on M Q nor 
on uj . 

It should be emphasised that in this subsection and hereafter the symbols a , e , uj , Q , i 
M , ft stand not for the exact values but for the mean values of the appropriate variables. 
A mean value of an element is understood to include the secular and long-period parts, 
the short-period components being averaged out. 

The case of uniform planetary precession ( p, = const ) was studied in Efroimsky 

(2005a). In that case, the terms containing fl evidently vanish. Besides, it turns out 
that, for constant /X , the mean values of the other /X-dependent terms, except one, van- 
ish too: 



df_ 

da 



x g - / x 



dg_ 



Cj = e , Q , u , i , M a , 



(26) 
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where 



A*± = A*i sin % sin f2 — /i 2 sin i cos f2 + yU 3 cos i 

(28) 

= Jp sin i sin f2 — /i p sin I p sin « cos Q + h p cos J p cos « 

is the projection of the planets' precession rate ft onto the instantaneous normal to the 
satellite's orbit. 4 

Hence, in this approximation and under the assumption of constant ft , in order to 
compute the secular parts of the orbital elements, it is sufficient to amend the Hamiltonian 
with the /z-dependent addition and to ignore all the other /X-dependent terms except 
the one given by (27). This will no longer be the case for variable precession, i.e., for 
time-dependent p, . Section 2 of our article will address itself to calculation of the secular 
parts (26 - 29) in the case of time-dependent p, . 



2 Equations for the first-order secular parts 
of the osculating elements. 

2.1 Two Fourier expansions of the precession spectrum 

Precession of the planetary spin axis has a continuous spectrum that spans from the polar 
wander and the fastest nutations to the Chandler wobble to the long-term variations whose 
time scales go all way to billions of years. When the planet has a sufficiently massive moon 
capable of influencing the planetary precession, the rate of this precession, /x , should be 
regarded as a function not only of time but also of the position of the satellite. We shall be 
interested, however, in the situation where the satellites are small and do not considerably 
influence rotation of their primary (while rotation variations of the primary still may affect 
the satellite orbits). This is, for example, the case of Mars whose tiny satellites affect its 
precession only in a very high order (Laskar 2004). Under these circumstances, it is fair 
to treat the precession rate as a function of time solely: 

fl(t) = r \p, {s \u) sin(wt) + ft (c) (u) cos(ut)l du , (29) 
Jo 1 J 

u standing for the angular frequency. In what follows, it will be convenient to describe 
the evolution not in terms of time but via the true anomaly v of the satellite. For our 
present purposes, it will be advantageous to express the precession rate as function of the 
satellite's true anomaly: 

jt(v) = f°° \fl {s) (W) sm(Wu) + fL (c) {W) cosiWv)] dW . (30) 
Jo 1 J 



4 Here /i± is expressed in the basis x , y , z associated with the planet's equator of date. Unit vector 
z is perpendicular to the equator of date, while x is pointing along the line of the ascending node of the 
equator of date on the equator of epoch; therefore, the components /i 3 are given by (3). 
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W being the circular "frequency" related to the true anomaly v . Evidently, p,(t) fZ(u) , jl(u) , 
and pi{W) are four different functions. We nevertheless denote them with the same no- 
tation p,( . . . ) because the argument will always single out which particular function we 
mean. The interconnection between functions f2(v) and jl{t) is given by 

ji(t) = fl(u)\ . 

is — | n at 

The interconnection between the Fourier components is less obvious. However, it simplifies 
under the assumption of vanishing eccentricity and slowly- changing semimajor axis: 

jt(W) « nft(u) | w = u/n , n = (Gm^a-^ 2 . (31) 
A rigorous relation to be used below is 5 : 

dft(u) _ dft(t) (dt\ _ (dt_^ (dM^ = /?(l-e 2 ) 3/2 2 ^ (32) 



dv dt \ d »Ja,e,i,u,,n,M \ dM J a ,... V ih ' n(l + ecosi/)' 



2.2 The role of the (p, •(...)) and (p, •(...)) terms 

These terms, ignored in the literature hitherto, implement the subtle influence of the 
planet's orbit precession upon its satellites' motion. The physical content of this effect is 
as follows: first, the precession of the planetary orbit slowly alters the solar torque acting 
on the planet; second, the variations of this torque entail changes in the planetary spin 
axis' precession; and, finally, third: these changes influence the satellites' orbits. This 
three-step interaction is extremely weak; still, its effect may accumulate over very long 
periods of time. 



2.2.1 The (...)) terms 

To illustrate the role of commensurabilities between the satellite orbital motion and the 
planetary nutations, let us consider the average ( /7 ■ ( (df/de) x g — (dg/de) x / ) ) 
emerging in the equation (22) for duo/dt and in the equation (25) for dM Q /dt : As shown 
in section A. 4 of the Appendix to Efroimsky (2004), 

_ I df _ <9g \ na 2 (3e + 2 cosu + e 2 cos u) , , 

yde de J (1 + e cos v) Vl - e 2 

/i± being given by (29). By virtue of (75) and (30), its secular part at some v will be: 

df _ =• dg\ 1-e 2 2 r' =w i x3e + 2 cos(u + u') + e 2 cos(u + u') 



^ \de B J dej 1 2tt h<=-^^ (1 + e cos(i/ + 1/')) 



2tt 



n a 



y y di/' [//^W sin(W(i/ + !/')) + A*±(W) cos(W(i/ + !/'))] [ (2 + e 2 ) cos(z/ 



rf^ (l 



5 This relation immediately follows from the well known equality dM (1 + e cosi/) 

,2\ 3 / 2 



, one upon which also the averaging rule (75) is based. 



11 



/ 5 \ 5 

+ u') + f - 3 e - - e 3 J cos2(z/ + z/') + 3 e 2 cos3(v + v') - - e 3 cos4(z/ + z/') + 0(e 4 ) 



(2 + e 2 ) ^ c) (l) + (-3e- \ e 3 ) ^(2) + 3e 2 ^ c) (3) - ^^(4) + 0(e 4 ) 



W being the angular "frequency" related to the true anomaly v , as in equation (30). Not 
surprisingly, the integral over W has been reduced to an infinite sum over the discrete 
values W = 1,2,3, 4, . . . corresponding to commensurabilities between the orbital 
frequency of the satellite and the nutational frequencies of the oblate planet. 6 The main 
resonant input comes from the principal commensurability W = 1 , i.e., from the nutation 
mode resonant with the orbit. The higher-order resonant inputs originate from the faster 
nutations characterised by W = 2 , 3, 4 . . . In equations (20 - 25), almost all terms 
proportional to ft produce such resonances. At the time when we are writing this paper, 
our knowledge about the fast nutations and polar wonder of Mars is yet very limited, and 
we shall not venture to offer quantitative estimates of the time scale over which the effect 
of these resonances upon the satellite orbit becomes considerable. 

Slower than W = 1 variations of fx bring no nonresonant contributions into the 
average of the right-hand side of (33). It can be shown that none of the ( fx- (...) ) terms 
emerging in (20 - 24) yield a nonresonant input. (Hence, (26).) For these reasons, in the 
rest of this paper, the terms (fx (•••)) will be omitted. 

2.2.2 The (fx •(...)) terms 

Let us consider, as an example, the term (fx ■ ( — / x (df/du)^j ) showing up in the 
right-hand sides of equations (21) and (23). We have from Appendix All of Efroimsky 
(2004): 

\ ou (1 + e cosz/) 

Just as in the preceding example (34), orbital averaging of this expression would yield res- 
onant terms entailed by commensurabilities between the orbital frequency of the satellite 
and the fast variations of fx . For the reasons explained above, here we omit these contri- 
butions. However, in distinction from the ( /x •(...) ) terms, some of the ( fx •(...) ) 
do have nonresonant components. For example, the mean part of (35) will be finite even 

for a constant fx : 

• / df \ , 2/ 2 \2(l-e 2 ) 3/2 r dv a 2 / 2 x 



6 It should be emphasised that (34) was obtained by a certain approximation: the averaging ignored 
the back-reaction of the short-period motions upon the long-period ones (i.e., it ignored the fact that, 
after each orbital period, the satellite does not return to exactly the same point it started); for example, 
it was assumed that the elements e and a remained constant during the integration over v' from 
through 2ir . 
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the superscript dot denoting a time derivative taken in the frame co-precessing with the 
equator of date. In other words, fi± is, by definition, not a full time derivative but its 
projection onto the instantaneous normal to the satellite's orbit. So defined fi± contains 
only derivatives of jij but not of the angles: 

= Ai sin i sin Q — fi 2 sin % cos Q + {13 cos i . (37) 

As shown in the Appendix below, fi± can be expressed via the longitude of the node, 
hp , and the inclination, I p , of the equator of date relative to the one of epoch: 

fi± = Ip sin i sin f2 — ( h p sin I p + h p I p cos I p J sin % cos f2 + ( h p cos I p — h p I p sin I p J cos 2 



~ /ip ( — sin Ip sin i cos Q + cos J p cos i ) . (38) 

The quantities h p , J p and their time derivatives can be calculated from integration of 
the Colombo equation of spin precession in inertial space, 

a (k- n) (k x n) , (39) 



dt 



T 



k = ( sin Ip sin /ip , —smI p cosh p , cos/p) being a unit vector pointing along the 

major-inertia axis of the planet, and n = ( sin I orb sin fl orb , — sin/ orfe cosQ or b , cos/ or b) 
being a unit normal to the planetary orbit plane defined (relative to some fiducial plane) 
through the inclination I or b and longitude of the node Q or b . The constant (or, better to 
say, the slowly-varying factor) a is given by 

2 3 (l - 4) c 

where n p , e p , s , and A , _B , C are the mean motion, the eccentricity, the spin angular 
velocity, and the moments of inertia of the planet. (As ever, C > B > A .) Even in the 
relatively simple case of C > B = A , the planet's axis of rotation does not describe a 
circular cone because the unit normal to the planet's orbit, n , is subject to variations 
caused by the precession of the planet's orbit about the Sun. While integration of the 
Colombo equation is explained below in the Appendix, here we would emphasise that this 
equation describes the evolution of planetary spin only under a very strong assumption of 
this spin being principal, i.e., in neglect of the Chandler wobble and polar wander. 



3 Evolution of the elements in the leading order of e 
3.1 The semimajor axis and the eccentricity 

As explained in subsection 2.2.1, the ( \x •(•••) ) terms may be omitted. The expressions 
for the orbital averages of the /x-dependent terms, derived in the Appendix, have the form: 

( h ■ ( - f x Mr ) ) = - ^ a 2 , (41) 



13 



-fx 



d_l 

duj 



— fi± a 



(42) 



Here 

fj,± = p, ■ w = Hi sin i sin Q — ji 2 sin i cos Q + /i 3 cos i , (43) 

where the unit vector 

w = x sin i sin Q — y sin z cos f2 + z cos i 

is the unit normal to the instantaneous plane of orbit, while the unit vectors x, y, z 
denote the basis of the co-precessing coordinate system x , y , z . (The axes x and y 
belong to the planet's equatorial plane, and the longitude of the node, Q , is measured 
from x .) 

The quantity fi± is defined as 



[l ■ w 



(44) 



but not as d(fl ■ w)/dt - a subtlety important to our further developments. 
Insertion of these expressions into equations (20 - 21) will give: 



da 



. = - 2a — VT 

dt n 



e 2 = 



de 
di 



5 fi± 



vT 



e 2 = 



2 a 5 / 2 
5 a 3 / 2 e 



vT 



2 n " 2 v/G^ 

For small eccentricities, the approximate solution is: 



(45) 
(46) 



a = a Q exp 



-2/3 



+ 0(e 2 ) 



1 + 



3 n n 



(47) 



e = e D exp 



-5/4 



+ 0(e 2 



1 + 



2 n„ 



,(48) 



where n G = (Gm) 1 ' 2 a^ 3 / 2 . We see that variations in the primary's precession exert 
almost no influence upon the satellite's semimajor axis and eccentricity. 

It should, nevertheless, be kept in mind that the satellite orbital elements evolve not 
only under the influence of the primary's precession but also under the action of tides. 
Within the truncated model developed in this paper, we shall neglect the tides, but shall 
introduce them on a subsequent stage of the project. 



3.2 The periapse, the inclination, and the node — in the leading 
order of e . 

Under the assumption of a and e remaining virtually unchanged, equations (22 - 24) 
will make a closed system, provided we omit the (p, ■ (...)) (for the reasons explained 
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above) and also substitute the orbital averages of the /x-dependent terms with their 
approximations in the leading order of the eccentricity. This level of approximation would 
be consistent with the approximation used in (47 - 48). As shown in the Appendix, 



fx 9 -l 

f X de 



, 



(49) 



fx ^ 

du 



= — a 2 ( fii sin i sinfi — p 2 sin i cosfz + p 3 cosi ) + 0(e 2 ) , (50) 



a 
~2 



-pi sin i sin fz cos i + p 2 sin i cos fz cos % — p 3 (2 — sin 2 ij + 0(e 2 ) , (51) 



di 



a 2 

= ( fii cosQ + fi 2 sinf2 ) + 0(e 2 ) 

2 



(52) 



Substitution of (19) and of the above expressions for the /i-terms into (22 - 24) will give 



us: 



dui 3 n J 2 (p. 
~dt = 4 



—) (5 cos 2 i - l) + u n cotz - u ± + - ( ^ cosfi + ^ sinfi) + 0(e 2 ) , (53) 
a/ v ^ 2\n n J 



di . fl± . p n 1 p 3 , 2 . 

— = — pi cosiZ — p 2 smiZ — — cot i — 1 + U(e ) , 

at n 2n sin i n 



(54) 



dVt 3 /p e 

— — = n Jo — 

dt 2 U 



Pn 

cost — - — : + 



sin i 2 sin i 



1 Pi n . P2 . n 

— — cosiZ H smiZ 



n 



where p± and are given by (43 - 44). The quantity 

p n = — pi sin cos i + p 2 cos fz cos z + p 3 sin z 



+ 0{e 2 ) 



(56) 



is the component of p, , pointing from the gravitating centre towards the ascending node 
of the orbit, while 



Pr, 



— pi sin Q cos % + p 2 cos Q cos i + p 3 sin i 



(57) 



is its time derivative taken in the frame co-precessing with the satellite orbit plane. (Taking 
the derivative in this frame, we differentiate only the components of p,, but not the angles.) 

Under the assumption of constant a and small e , equations (54 - 55) make a closed 
system. 



(55) 
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3.3 Goldreich's approximation 



It would now be tempting to introduce an even stronger assumption that both | ft 
|/ ( n 2 J 2 sin i ) and \f2\/ (nJ 2 sin i ) are much less than unity, and to derive therefrom 
the system 

rfr 7 ~ " 2 nJ2 UJ (T^if ' (58) 



di 

whose solution, 



— Hi cosf2 — f/,2 sinf2 , (59) 



i = - — cos [ - x (t - Q + Q ] + — sin [ - \ (t - t Q ) + Q G ] + i , 

_ , \ ^ , 3 T /p e \ 2 cosi 

" = — y (t — t ) + u where y = - n Jo — o , 

AV ° ; A 2 2 Vay(l- e 2) 2 ' 



(60) 



seems to indicate that, in the course of planet precession (the term "precession" including, 
as agreed above, also nutations and the Chandler wobble and polar wander), the satellite 
inclination oscillates about i . Approximation (60) has already appeared in the literature. 
Goldreich (1965) derived such equations for the orbital averages of some nonosculating 
elements (which later were termed, by Brumberg (1992), "contact elements"). In our 
case, however, the approximation (60) was derived for the secular parts of osculating 
elements. We see that, in neglect of /x 2 -terms and under the assumption of constant 
p, , the equations for the secular parts of osculating elements coincide with those for the 
secular parts of the contact ones. (For a detailed explanation of this fact see Efroimsky 
(2005a).) 

The evident flaw of approximation (58 - 60) is its invalidity in the closest vicinity of 

the equator. In this vicinity, the parameters | ft |/ ( n 2 J 2 sin i ) and \p,\/ (n J 2 sin i ) are 
no longer small; so the entire approximation falls apart and gives no immediate indication 
on whether the inclination will go through zero and alter its sign or will "bounce off" the 
equator. At the first glance, this technical subtlety does not affect this approximation's 
main physical outcome, one that the inclination remains limited and shows no secular 
increase. In reality, though, the matter needs further exploration. For example, if the 
orbit keeps bouncing off the equator and the sign of i stays unaltered for long, then 
the term jj, n / (2n) in (54) may, potentially, keep accumulating through aeons, creating a 
drift of the inclination. Whether this is so or not can be learned numerically through a 
more accurate approximation based on equations (53 - 55). A more definite thing is that 
the Goldreich approximation is intended only for low inclinations: as can be seen from 
equation (55), at high inclinations it will fail, because the term fij_/ sin i will dominate 
over the J 2 cos i term. All these issues will be addressed in our subsequent paper (Gurfil, 
Lainey & Efroimsky 2007). 

3.4 Can precession cause secular changes of the inclination? 

Above we saw that the Goldreich approximation reveals no secular terms in the expression 
for the inclination relative to the moving equator. While a reliable quest into this matter 
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will demand numerical integration of the entire system of the planetary equations, we shall 
try to work out a qualitative estimate based on the approximation less crude than that 
of Goldreich. To this end we shall plug (44) and (56) into (54), and shall omit all the 
long-period terms. Thus we shall be left with the following estimate for the secular part 
of i : 

, .(sec) 

CLl 1U3 . o\ / \ 

— - — = — sin! + U(e ) + long-period terms (61) 
dt 2n 

3.4.1 Small initial inclinations 

For small inclinations, the above equation will look: 



, .(sec) 

di 
dt 

where 



A i i3ec) (62) 



A _ £3 _ h p cos/ p _ 

2n 2n 2n ' K J 

h p and I p being the longitude of the node and the inclination of the equator of date on 
that of epoch, (see Appendix A8). From here we see that the osculating component of i 
will, approximately, obey 

t sec> « i e At . (64) 

The exponential dependence evidences of the presence of chaos in the system. It should 
be mentioned, though that the chaos will be weak, because A is extremely small. Besides, 
the second derivative of the precessing equator's node, h p , which enters the expression 
for A , does not keep the same sign through aeons. The rate, at which node h p and 
the inclination I p evolve, can be computed, via the Colombo equation, from the rate of 
precession of the planet's orbit about the Sun. (For details on the Colombo model see the 
Appendix below.) Qualitatively, one may expect the spectrum of h p and I p to resemble 
the frequency content of the planet's orbital precession. (See the table in the Appendix.) 
On all these grounds, the time dependence of i is constituted by the high-frequency 
oscillations (60) superimposed on a much slower evolution (64). We expect this slow 
evolution to look as a saw-tooth plot, because the sign of h p (and therefore of A ) alters 
from time to time for the reason explained above. Due to this saw-tooth nature of the 
long-term evolution of i^ sec ^ , no considerable secular increase of the satellite's inclination 
should be expected, at least in the case of a small initial i Q . Numerical calculations 
performed in the e 3 order confirm this conclusion. Moreover, it turns out that even at 
not too small initial inclinations no secular changes in i accumulate over time scales of 
order billion years. 

The said numerical results and plots are presented in our subsequent paper Gurful, 
Lainey & Efroimsky (2007), devoted to a numerical implementation of our semianalytical 
model in the e 3 order. 



17 



3.4.2 Large initial inclinations 

For near-polar orbits, the equation (61) will read: 

. (sec) 



di 

— « A , (65) 

whence 

i sec) w At . (66) 

Once again, due to the undulatory sign alterations of A , we shall get a "saw-tooth" 
behaviour, though this time the teeth will be less steep than in the small-inclination case 
governed by the exponent (64). The teeth will be expected to cross the polar orbit once 
in a while. This kind of time dependence (so-called "crankshaft") is indeed what results 
from the numerical computations. (Ibid.) 

All in all, unless we begin very close to the pole, the variable equinoctial precession is 
not expected to entail secular changes in the satellite inclination relative to the equator of 
date. Exceptional is the case of near-polar orbits: in that case, leaps across the pole are 
possible. (See Ibid, for details and plots.) 

4 Preparation for computation in the e 3 order 

Insertion of (19) into equations 7 (20 - 24) will lead us to the following system: 
da fi± / 2 \i/2 

-=-2-a (1-e) , (67) 

de 5 li± / o \ 1/2 

e 



dt 2 n 



du 3 n.J 2 fp e \ 2 /5 2 . 1\ cos? • / - df , 

■ cos i - - -//± + A* n cotz — , — : (/x -/ x — ) , (69) 



dt 2(l-e 2 f\a) \2 2 J ^ ^ n na 2 (l - e 2 ) 1 / 2 sin % \ J dt 



di 
dt 



Hi cosf2 — /i 2 sinfi 



cosi . • / r df \ . 1 . • ( ? df . 

+ ^71 2M/2 • ■ ( /"* ~/ X 77~ ) 271 2U/2 • ■ ^ \ ~f x 7^7 ' 70 

na 2 (1 — e 2 ) 1 / 2 sin « \ aa; / na 2 (1 — e 2 ) 1 / 2 sin % \ <7i2 ' 



7 As within this model both the Hamiltonian perturbation and the /2-dependent terms are substituted 
with their orbital averages, M a becomes a nuisance parameter, so the planetary equation for dM D /dt is 
omitted. 
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dQ 3 

~dt ~ ~2 U 2 I a J (1- e 2) 2 sinz ' n a 2 (1 - e 2 ) 1 ^ s in % 



+ 



where, according to the Appendix, 
<*■[-/ xf 



a 

T 



| /ix — (2 + 3e 2 ) cosf2 + 5e 2 (cosf2 cos2c<j — sinf2 sin2cj cosi) 



+ 



A*2 



— (2 + 3e 2 ) sinf2 + 5e 2 (sinfi cos2cj + cosf2 sin2cj cos i) 



+ 



A*3 



5 e 2 sin 2w sin z J j 



(72) 



a 2 

— (^2 + 3 e 2 ) ( fii sin i sin f2 — /i 2 sin % cos f2 + /i 3 cos i ) , 



Of. 

— I /ii sm z 



(2 + 3 e 2 ) sin f2 cos i + 5 e 2 ( cos f2 sin 2c<j + sin f2 cos 2c<j cos i 



+ /i 2 sin % (2 + 3e 2 ) cosf2 cos i + 5e 2 (sinQ sin 2a; — cosf2 cos2c<j cos i) 



(2 + 3 e 2 ) (2 - sin 2 i ) + 5 e 2 sin 2 i cos 2cj ] } 



5 Conclusions 



In this article we continued our analytical investigation of the behaviour of orbits about 
a precessing oblate planet. We built up a reasonably simplified model that takes into 
account both the long-term variability of the planetary precession (variability caused by 
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the planet's orbit precession) and the short-term variability (polar wonder, etc). 

We have written down equations (67 - 71) that describe evolution of the satellite orbit at 
long time scales. These equations also include known functions of time, ^ , /i 2 , Hz , /x± , \i n , , 
fi 2 , A3 , ft± , f^n , which are various projections of the planet axis' precession rate. An 
algorithm for computation of these functions of time is presented in the Appendix. These 
functions vary in time as a result of precession of the primary's orbit about the Sun. This 
way, we have analytically established connection between the precession of the planet's or- 
bit and the evolution of its satellites. Physically, this connection comes into being through 
the following concatenation of circumstances: precession of the planetary orbit leads to 
variations in the Solar torque acting on the planet; the torque variations cause changes in 
the planet axis' precession; these changes, in their turn, entail variations of orbits of the 
planet's satellite. This effect is extremely weak and accumulates over very long time scales. 

All in all, we have fully prepared a launching pad for computation of the evolution of 
near-equatorial circummartian orbits at long time scales. The methods and results of this 
integration, and their physical interpretation will be presented in our next publication 
(Gurfil, Lainey, & Efroimsky 2007), which will also include the pull of the Sun. Briefly 
speaking, those results are to be three- fold. First, it turns out that our semianalytical 
model is robust beyond expectations. Despite the averaging and the neglect of /z 2 -terms, 
it works very well over timescales up to, at least, 20 Myr. Second, it turns out that 
precession by itself (i.e., in the presence of the Sun but in the absence of the other phys- 
ical factors like the tides and the planet's triaxiality) cannot cause accumulating secular 
changes in the satellite inclination, provided the initial inclination is not too large. This 
means that, for orbits not too close to the polar one, the main prediction of the Goldreich 
model stays valid, even though the model cannot adequately describe the entire dynamics 
(which becomes weakly chaotic). Third, it turns out that in the vicinity of the polar or- 
bit precession of the primary can cause major alteration of the satellite orbits, including 
unusual features in the behaviour of the inclination. See Ibid for more details. Further 
work along this line of research will be aimed at including more factors into the model - 
the tidal forces, the pull of the Sun, the triaxiality of the planet, etc. 
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Appendix. 

The goal of this Appendix is to calculate, in neglect of nut at ion- caused resonances, the 
secular parts of the /T-dependent terms emerging in the planetary equations (19) - (25). 
We shall also explain how to compute the time dependence of various projections of the 
planetary precession rate ft and of its time derivatives. 
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A.l. The averaging rule 



The mean values are to be calculated via the averaging rule: 

2 TT J-n (1 + e COSU) 

Since the averaging is carried out over the true anomaly, it will be convenient to express 
the precession rate not as a function of time, fX(t) , but as a function of the true anomaly: 

jt(v) = r \ft (s) (W) sm(Wu) + ft {c) (W) cosCWv)] dW . (76) 
Jo 1 J 

W being the circular "frequency" related to the true anomaly v . 

In what follows, we shall need the average of the projection of fi{t) onto the instan- 
taneous normal to the orbit. This projection, fj,± , will be expressed by 

fj,± = p, ■ w = fj,i sin i sin Q — /i 2 sin % cosf2 + /i 3 cos i . (77) 

where the unit vector 

w = x sin % sin Q — y sin i cosQ + z cos i (78) 

stands for the normal to the instantaneous osculating ellipse, and the unit vectors x , y , z 
are the basis of the co-precessing coordinate system x , y , z . (The axes x and y belong 
to the planet's equatorial plane, and the longitude of the node, Q , is measured from x .) 
Expressions of fij via the longitude of the node and inclination of the equator of date 
relative to that of epoch are given below in section A. 12. 



A. 2. Calculation of the secular part of u ■ I - f x ^ 

y ' oa 

According to the Appendix to Efroimsky (2004), the term of our concern, written in 
terms of the orbital elements, is given by 

i* ( ~ ^ X fa") = \ ^ a U ^ 1 ~~ to ^ ^ l _ ^ ' 

expression linear in t — t Q . It coincides with its secular part, if we assume that the 
spectrum of ft lacks short-period modes. 



— df 

A. 3. Calculation of the secular part of jt • I - / x 



de 



From Efroimsky (2004) we have: 



H 1 -/ x 9 -f ] = -fi ± a 2 -g ^ sinz, . . (79) 
oe I 1 + e cos z/ 
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The secular part of this expression, calculated through 

.2^3/2 

' / - 

+ e cos v) 



, '-> ( ^ df\ , 9 /, 9 \ (1 — e 2 ) 3//2 f n sinz/ dv . ^ N 

{it ■ [-fx -£)) = -H±a 2 1-e 2 V J / — -3 , .(80) 

\ oe v ' z 7r J-7T 1 + ecosz/ 



will vanish, because the function in the numerator is odd, while the expression in the 
denominator is even. 



A. 4 Calculation of the secular and long-period parts of ft ■ y - f 

We have from Efroimsky (2004): 

jt. (-fx 9 /) = - > ± a' {1 ~ ^ . (81) 

y au J (1 + e cosz/) 

This expression has the following secular part: 

< h ■ ( - f x 3 i) > = - Ax a* (1 - e f <i^2!! r *> = , 82) 

\ aa; / v ' 2 tt ./-tt (1 + e cosz/) 



. 2 (1 - e 2 ) 7/2 2 + 3 e 2 a 2 / 2 x 

—^h^jr^r = -^y ( 2 + 3e ) = 

a 2 

— (2 + 3 e 2 j ( fii sinz sinf2 — fi 2 sin z cosf2 + /i 3 cosz) 

Be mindful that the time derivative of fi± , calculated with aid of (77), reads as 

fi± = A*i sin i sinQ — /t 2 sin z cosf2 + /t 3 cosz (83) 

and contains only derivatives of /i-,- but not of the angles. This happens because fi± is 
not a full time derivative but is defined through fi± = p, ■ w . 

A. 5. Calculation of the secular and long-period parts of p, ■ ^ - f x 

In this subsection and below we shall need the following auxiliary integrals: 
r [ n 1 a 2 + 3e 2 

T ° = / TTL ^ dv = 7v , (84) 

J-n (1 + e cosz/) (1 — e 2 ) 1 

„ /•*" sin(w + z/) cos(^ + v) 5 2 sin(2^) 

Tl = / 7T1 ^ ciz/ = - tt e , (85) 

J-tt (1 + e cosz/) 2 (1 — e 2y/ 
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T 2 = / — dv 



7T (1 + e cosu) 



1 2 + 3e 2 - 5e 2 cos(2o;) 
2^ (1 - e 2 ) 7 / 2 



(86) 



The first component: 



(87) 



2N 2 



( /ti a 2 



(1 - e 2 ) 



(1 + e cosz/)' 



[ cos f2 cos(a; + v) — sin Q sin(a; + z/) cos i ] sin(a; + z/) sin z 



■ 2 2A 2 (1 - e 2 ) 3/2 

= ^ a (1 - e ) 2?r 



{ Ti cosf2 sin i — Y 2 sinQ cos z sin z } 



a 2 

/ii — sin z | — (2 + 3e 2 ) sinf2 cos z + 5 e 2 [cosf2 sin(2cj) + sinf2 cos(2cj) cosz] | 



The second component: 



( /i 2 a 2 



2^ 



1 - e 2 ) 



(1 + e cosz/)' 



[sinfi cos(u; + v) + cosf2 sin(a; + v) cos z ] sin(u; + z/) sin z 



/i 2 a 2 (l - e 2 ) 



2 (1 - e 



2\3/2 



2tt 



{ Ti sinfi sin z + T 2 cosf2 sin z cos z } 



/i 2 a 



3 vl ~ e 2 ) 7/2 
2tt 



sin (2a;) 
2 " (1 - e 2 ) 



- 7r e 



2^/2 



sinf2 + 



1 2 2 + 3 e 2 - 5 e 2 cos(2a; 



- 7r e 
2 



(1 - e 2 ) 7/2 



cos Q > sin z 



a 2 

= /i 2 — sin z | (2 + 3e 2 ) cosf2 cos z + 5e 2 [sinfi sin(2a;) — cosf2 cos(2a;) cosz] j 



The third component: 



= ( - A*3 « 2 



2^ 



(1-e, 
'1 + e cos i/) 



cos 2 (a; + z/) + sin 2 (a; + z/) cos 2 z ) = (89) 
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-/i 3 a 



:i-e 2 ) 

2 7T 



2N7/2 



(t -T 



2 sm 2 i) = -/i 3 a 2 



1 - e 2 ) 7/2 f 2 + 3e 2 2 + 3e 2 - 5e 2 cos(2a; 

< 7T =777 — 7T ^ 



2 7T 1 (l-e^) 7 ^ 



2 (1 - e 2 ) 7/2 



sin i 



2 — sin i 



+ 5e 2 sin 2 i cos(2cj) | 



Total: 



I fix sin i — (2 + 3 e 2 ) sin f2 cos z + 5 e 2 ( cos f2 sin 2u + sin f2 cos 2uj cos i ) 



+ fi 2 sin i (2 + 3e 2 ) cosf2 cos z + 5e 2 (sinf2 sin2cj — cosf2 cos2c<j cos i) 



~ fJ>3 



(2 + 3e 2 ) (2 - sin 2 z) + 5e 2 sin 2 i cos2cj] } 



(90) 



Interestingly, even in the limit of vanishing eccentricity this sum survives and becomes 



— j — fii sin i sinf2 cos i + /i 2 sin i cosQ cos z — /i 3 (2 — sin 2 | = — fi n sin i — a 2 /i 3 . (91) 



Moreover, even when both the eccentricity and inclination are nil, this sum still remains 
nonzero and is equal to — /i 3 a 2 . 



A. 6. Calculation of the secular and long-period parts of fx ■ I - / x 

<)i 



The first term: 



df 

di 



(92) 



( /ii a 2 



2^2 



(1 - e 2 ) 
[1 + e cosu) 



[ — cosfl sin(cj + v) — sinf2 cos(cj + u) cos i ] sin(c<j + v) ) — 
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2 fl 6 2 )^^ 

a 2 (l — e 2 ) [ — T 2 cosf2 — T x sinf2 cos i 



(l-e 2 ) 7/2 



2tt 



vr 2 + 3e 2 - 5e 2 cos(2cj) 
2 (1 - e^) 7 / 2 



cosf2 — 



5 9 sin(2u;) 

— 7T P — 

2 (1 - e^) 7/2 J 



sin f2 cos i 



a 2 

= /ii — | — ^2 + 3e 2 ) cosf2 + 5e 2 [cosf2 cos(2c<j) — sinfi sin(2cj) cos i 



The second term: 



(93) 

(1 - e 2 ) 2 

/i 2 a 2 2 [ — sinf2 sin(a; + z/) + cosf2 cos(cj + z/) cos i ] sin(cj + v) ) — 

(1 + e cosu) 



fi2 a 2 (l - e 2 ) 



2 (1 - e 2 ) 



2\3/2 



27T 



[ — T 2 sin + Ti cos Q cos z ] = 



A*2 a 



(1-e 



2\7/2 



27T 



vr 2 + 3e 2 - 5e 2 cos(2cj) 
2 (1 - e^) 7 / 2 



sinfi + 



5 2 sin(2o;) 

- 7r e =j 

2 (1 - e 2 ) 7/2 . 



cos Q COS z 



a 2 

— A2 — { — (2 + 3e 2 ) sinf2 + 5e 2 [sinfi cos(2cj) + cosf2 sin(2w) cos? 



The third term: 



(A, if */ 



,2\2 



( A3 a 2 



_n 

(1 + e cos z/) z 



sin(cj + z/) cos(u; + v) sin i 
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Ms 



2 (l-e 2 ) 7/2 5 2 sin(2cj) 

- , sin i = u 3 a - 7r e =7^ sin z 

2 7T ^ 2 7T 2 (1 _ e 2) 7 /2 



/i 3 a 2 - e 2 sin 2a; sin i 



Total: 



(95) 



cr r 



— (2 + 3 e 2 ) cos Q + 5 e 2 (cos f2 cos 2u; — sin Q sin 2u; cos z) 



+ 



1^2 



(2 + 3e 2 ) sinfi + 5e 2 (sinfi cos2cj + cosf2 sin2cj cosi) 



+ 



A*3 



5 e sin 2a; sin i 



In the limit of vanishing eccentricity this sum approaches 



77 { Ai cos ^ + A2 sinfi } , 

expression that bears no dependence upon the inclination. 



(96) 



A. 7. Calculation of the secular part of ft • - / x 



df 



The considered quantity, when expressed through the orbital elements, looks as 

df 



M • - / x 



— fi± a 2 y/l — e 2 



(97) 



As we are trying to calculate not the exact values of the elements but their secular parts, we 
assume that, at shorter than an orbital period time scales, a , e and fi± stay unchanged. 
In this approximation, the above expression coincides with its secular part: 



<P ■ [ ~f x 



df 



— fi± a 2 y/l — e 2 



(98) 
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A. 8. The planetary precession rate fl and its projection fi± 
onto the satellite's orbital momentum. 



Let the inertial axes ( X , Y , Z) and the corresponding unit vectors ( X , Y , Z ) 
be fixed in space so that X and Y belong to the equator of epoch. A rotation within the 
equator-of-epoch plane by longitude h p , from axis X , will define the line of nodes, x . 
A rotation about this line by an inclination angle I p will give us the planetary equator of 
date. The line of nodes x , along with axis y naturally chosen within the equator-of-date 
plane, and with axis z orthogonal to this plane, will constitute the precessing coordinate 
system, with the appropriate basis denoted by ( x , y , z ) . 

In the inertial basis ( X , Y , Z ) , the direction to the North Pole of date is given by 



sin I p sin h p , — sin I p cos h p , cos I p ) 



(99) 



while the total angular velocity reads: 

->(inertiaX) 



total 



= Z S + ^ inertial ^ 



(100) 



the first term denoting the rotation about the precessing axis z , the second term being 
the precession rate of z relative to the inertial frame ( X , Y , Z ) , and s standing for 
the angular velocity of rotation about the axis z . This precession rate is given by 



// "' = ( cos hp 



-*(inertial) 



I p sin h p , h p ) 



(101) 



because this expression satisfies ^ %nertiaV ) x z = z . 

In a frame co-precessing with the equator of date, the precession rate will be represented 
by vector 



R 



(inertial) 



(102) 



where the matrix of rotation from the equator of epoch to that of date (i.e., from the 
inertial frame to the precessing one) is given by 



— >p 



cos h r 



cos Ip sin/ip 



sin Ip sin h p 



sin h r 



cos I p cos hp 



sin Ip cos hp 







sin I n 



cos Z„ 



(103) 



From here we get the components of the precession rate, as seen in the co-precessing 
coordinate frame (x , y , z) : 



H = { jii , /i 2 , /i 3 



(4 , 



h„ sin/„ 



hp cos I p 



(104) 
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In our paper we also need the components of ft dot standing for derivatives calculated 
in the frame co-precessing with the equator: 



/ii , /i 2 , A*3 ) = ( i P , hp sin I p + h p I p cos I p , h p cos I p - h p I p sin I p 



(105) 



The matrix of rotation from the precessing frame of the equator of date to the frame 
associated with the satellite's orbital plane will look: 



cosf2 



sinfi 







— cos % sin f2 cos i cos f2 sin i 

sin i sin Q — sin % cos fl cos i 

This will give us the precession rate as seen in the instantaneous orbit frame: 

/T (or6) = R p ^ it . 



(106) 



(107) 



The second component of this vector (i.e., the component pointing toward the ascending 
node of the satellite orbit relative to the equator of date) is what we need in our formulae 
(52) and (54): 



— — A*i sin ^ cos i + A*2 c °s ^ cos % + /x 3 sin i 



= — i p sin f2 cos z + h p sin 7 p cos f2 cos i + h p cos J p sin i . 



(108) 



Its time derivative (taken in the frame of reference co-precessing with the equator of date) 
is: 



— fii sin Q cos % + fi2 cos f2 cos i + {13 sin i 



= — i p sin f2 cos i + (^h p sin J p + /i p i p cos J p ) cos Vt cos z + ( h p cos J p — /i p i p sin J p ) sin i . 

The third component of fi!-"^ (i.e., the component orthogonal to the instantaneous plane 
of the orbit) is exactly what we need in (27) and (33): 



(109) 



A*± = A*i sin 2 sin f2 — /i 2 sin i cos Q + fi 3 cos i 



Ip sin i sin f2 — /i p sin I p sin z cos Q + h p cos J p cos i 



(110) 



Its time derivative fi± defined in the axes co-precessing with the equator (and therefore 
equal to ft ■ w , not to d(p, • w)/dt) will now be expressed by 

= fii sin i sin Q — fi 2 sin i cos Q + fi 3 cos i = 
i p sin « sinf2 — (Ji p sin/ p + /i p / p cos/ p ) sin i cosf2 + (/i p cos/ p — h p I p sin/ p ) cos i (HI) 



/i p ( — sin I p sin i cos Q + cos J p cos % ) , 



(112) 
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Ip and h p being the inclination and the longitude of the node of the equator of date 
relative to the one of epoch. 

The expression for fi± permitted approximations shown above because, for Mars' 
equator, the speed of the nodes' motion, \h p \ « 360°/ (1.75 x 10 5 yr) xs 2 x 10~ 3 yr^ 1 , 
much exceeds the rate of its inclination change, \I p \ w 5°/(0.5 x 10 e yr) rs 10~ 5 yr^ 1 
(Ward 1974). 



A. 9. Calculation of h„ and /, 



The question now becomes as to how to calculate the time dependence of h p and I p . 
As very well known, these two angles evolve in time because a non-spherical planet behaves 
itself as an unsupported top whose precession is instigated by the solar torque. The torques 
produced by the satellites are irrelevant (Laskar 2004), the cases of the Moon and Charon 
being exceptional. When the moments of inertia of the planet relate as C > B = A , 
the solar torque is 

*= ^ (O-A) (f-k) (fxk) , (113) 
while in the general case of C > B > A it is equal to 

provided that the spin mode is not too deviant from the principal one, and that this spin 
is much faster than the planet's orbital revolution about the Sun. 

In the above two formulae, M is the solar mass, R denotes the distance between 
the centres of masses of the planet and the Sun, the unit vector f points from the planet 
toward the Sun, and the unit vector k points in the direction of the major-inertia axis 
of the planet. 

The precession of the angular momentum L of the planetary spin obeys 

Colombo (1966) averaged this equation over the planet's year, under the assumption that 
the perturbing torque causes only very small variations of spin. This averaging yields the 
Colombo equation valid at timescales much exceeding one year: 



(! - 4)' 



n p and e p being the mean motion and the eccentricity of the planet's orbit about the 

— * 

Sun, and n being the unit vector normal to the planetary orbit; while dL/dt should be 
understood as a change of L over a year, divided by the length of the year: AL / P . The 
angular velocity of the planet about its axis being denoted with letter s , the Colombo 
equation may be rewritten as 
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the factor a being defined as 



a = 



3n 2 



C - 



A + B 



2 s (l - e£) 



3/2 



c 



2tt D 



(» 



o2A 



3/2 



c - 



A + B 



C 



(118) 



where P = 2njn v is the duration of the planet's year, and D = 2ix/s is that of its 
day. The relative difference between the moments of inertia may be expressed through the 
parameter J 2 emerging in the expression for potential via the planetocentric latitude 



V = 



Gm 



1 " E J n (jj ^n(sin0) 



[ Pe\ n 

2 S I ~ ) ^nj(sin0) cosj(A - A„j) ,(119) 

n=2 j=l V r y 



(where p e stands for the mean equatorial radius of the planet): 



Jo 



C - 



A + B 
2 



C - 



A + B 
2 



C 



M pi C M p] 

It is also interconnected with the nonsphericity parameter J : 



(120) 



J = 




3 g - j-^j C 
2 C Mp 2 



3 g - 

2 C 



AT 



(121) 



where p is simply the mean (not the mean equatorial) radius of the planet, while the 
quantity 



K = 



C 



m p z 



(122) 



is the squared ratio of the gyration radius \JcJm of the planet to its mean radius p . 
We thus see that 



C - 



A + B 
2 



C 



whence 



2tt D 



J 



J_ 

K 

2tt D 



3 mpl 
2 J2 — 



a 



P P (1 - ef 2 K P P (1 - e 2 ) 3/2 2 

2.95 x 10~ 3 



C 



(123) 



(124) 



which gave him the 



Ward (1974) used, for Mars, K = 0.359 and J 
value: 8 a Mar3 = 1.26 x 10~ 12 rad/s = 8.19 arcsec/yr . 

To further simplify the above expression (117), Colombo (1966) assumed that the spin 

— * 

angular momentum L is parallel to the spin angular velocity s : 



Cs k C s k 



8 In his formulae, Ward (1973) missed or deliberately neglected the factor (l — 



-3/2 



(125) 
In the case 

of Mars, this may look legitimate because nowadays this factor amounts to 1.013 . However, the Martian 
eccentricity is wont to have varied through aeons within the interval of e = 0.01 -j- 0.14 (Murray et 
al 1973). This means that the said factor, (l — e^) , might have varied from almost unity through 
1.03 . This 3 % increase will look less than innocent if we recall that several authors (Ward 1974, Laskar 
& Robutcl 1993, Touma & Wisdom 1994) insist on the stochastic nature of Mars' obliquity variations. 
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While investigating the dynamics of the Moon at less than cosmic time scales, Colombo 
certainly could afford this approximation. Compare the latter with the exact expression 
for L through s and through the moments of inertia C > B > A : 

L = \ Sl A + js 2 B + ks 3 C =\ Sl (A - C) + ]s 2 (B - C) + Cs (p - k) + C sk ,(126) 

s = (isi + j s 2 + ks 3 ) s" 1 being the instantaneous direction of the angular velocity of 
the planet's spin. We see that Colombo's approximation stems from the frivolous assertion 
that the planet always remains in a principal spin state. Indeed, insofar as s coincides 
with k the components S\ and s 2 are nil, and (125) becomes exact. Under such an 
assertion, the equation for unit vector aimed in the direction of the major-inertia axis, k , 
assumes the form 9 

^ = a (n ■ k) (k x n) (127) 

the unit normal to the planet's orbit, n , being subject to variations described by the 
formulae and tables worked out by Brouwer & van Woerkom (1950). 10 The quantities 
hp , I p and their time derivatives can be calculated from integration of the above equation 
for k . To this end, let us recall that our unit vector k coincides with the afore discussed 
unit vector z (see formula (99) above). Therefore, in the frame of the equator of epoch 
(which we assume, for convenience, to coincide with the ecliptic of 1950), k and dk/dt 
will read: 

~ T 

k = ( sin I p sin h p , — sin/ p cosh p , cos/ p ) (128) 

dk / • • • • \ T 

— = [Ip cos Ip sin hp + h p sin/ p cosh p , — I p cosI p cosh p + h p sin/ p sin/i p , — I p sinJ p J , (129) 

while the components of n , 

T 

n = ( sin I or b sin Q or t> , — sin I or b cosfl or b , cos/ or f,) (130) 
may be expressed through the auxiliary variables 

q = sin I or b sinQ or b , p = sin I or b cosQ or b (131) 
whose evolution will be found from 

oo 

9 = E N i sin + 5 i) > ( 132 ) 
i=i 

oo 

p = N j cos + 6 i) ■ (133) 

3=1 

Under the assumption that the orbital elements are defined relative to the ecliptic plane of 
1950, Brouwer & van Woerkom (1950) calculated the values of the amplitudes, frequencies, 
and phases used in the above formulae. Below follow the triples of numbers: 



9 By basing his research on the approximation (127), Ward (1974) implicitly made the strong assump- 
tion of Mars always remaining in the principal spin state, polar wander and nutations and the Chandler 
wobble being neglected. By employing a Hamiltonian, that generates equation (127), Laskar & Robutel 
(1993) and Touma & Wisdom (1994), too, rested their study on the same assumption. 

10 Brouwer & van Woerkom (1950) chose the ecliptic of 1950 as the reference plane. Since our eventual 
goal is to simply estimate the range of variations of i and £1 over large time scales, we can accept, 
without loss of generality, that at some distant epoch the Martian equator coincided with that plane. 
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j 


AT 


c> 

S 3 

(arcsec/yr) 


6 J 

(deg) 


1 


0.0084889 


- 5.201537 


19.43255 


2 


0.0080958 


- 6.570802 


318.05685 


3 


0.0244823 


- 18.743586 


255.03057 


4 


0.0045254 


- 17.633305 


296.54103 


5 


0.0275703 


+ 0.000004 


107.10201 


6 


0.0028112 


- 25.733549 


127.36669 


7 


- 0.0017308 


- 2.902663 


315.06348 


8 


- 0.0012969 


- 0.677522 


202.29272 



Technically, the computation of time evolution of I p and h p can be implemented 
through a set of differential equations obtained by substitution of 

T 

h — ( q , — p , u) , q = sin/ orb sinf2 orfe , p = sin/ orf) cosf2 or {, , u = cos/ or i, , (134) 

k = ( Q , — P , U ) , Q = sin/p smh p , P = sml p cosh p , U = cos/ p , (135) 
into the Colombo equation (127). Here follow these equations: 

- a [ q(t) Q{t) + p(t) P(t) + u(t) U{t) ] [ - p{t) U(t) + u(t) P(t) ] , (136) 
« [ q(t) Q(t) + p(t) P(t) + u(t) U(t) ] [ u(t) Q(t) - q(t) U(t) ] , (137) 

- a [ q(t) Q(t) + p(t) P(t) + u(t) U(t) ] [ - q(t) P(t) + p(t) Q(t) ] , (138) 
where, at each time step, the following values of q(t) , p(t) , and u{t) are to be used: 

oo 

q(t) = N j sin fa + Sj) , (139) 

oo 

pit) = N 3 C0S + S j) . ( 14 °) 

J'=l 

«(0 = ± - g(t) 2 - P (ty . (i4i) 

The resulting values of Q , P , C/ , obtained through this integration, will, at each time 
step, give us the angles I p and h p via formulae that follow from (134 - 135): 

h p = arctan^ , I p = arccosf/ , (142) 
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dQjt) 

dt 
dP(t) 

dt 

dU{t) 
dt 



It is evident from (135) that the variables Q(t) , P(t) , U(t) obey the constraint 

Q(t) 2 + P(t) 2 + U(t) 2 = 1 ; (143) 

and therefore fulfilment of this constraint should be checked during integration. Deviation 
from it will indicate accumulation of errors. At each step, some attention will be needed 
also when the current value of u(t) is evaluated. (We mean the choice of sign in (141).) 

Finally, it should be emphasised that the development by Brouwer & Clemence (1950) 
is limited in terms of precision and, therefore, in terms of the time span over which it 
remains valid. A more accurate and comprehensive development, with a validity span of 
tens of millions of years, was recently offered by Laskar (1988). At the future stages of 
our project, when developing a detailed physical model of the satellite motion, we shall 
employ Laskar's results. 
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